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Abstract 



We first consider an exactly solvable classical field model to understand the coherence 
properties and the density fluctuations of a one-dimensional (ID) weakly interacting de- 
. generate Bose gas with repulsive interactions at temperatures larger than the chemical 

^ ' potential. In a second part, using a lattice model for the quantum field, we explain how 

q I to carefully generalize the usual Bogoliubov approach to study a degenerate and weakly 

interacting Bose gas in ID, 2D or 3D in the regime of weak density fluctuations. In the 
last part, using the mapping to an ideal Fermi gas in second quantized formalism, we 
^S) , calculate and discuss physically the density fluctuations and the coherence properties of 

^ ' a gas of impenetrable bosons in ID. 
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1 A CLASSICAL FIELD MODEL FOR THE ID WEAKLY INTERACTING BOSE GAS 

1.1 Reminder: quantum theory for the ideal Bose gas 

1.1.1 Second quantized formalism 



, We consider the case of spinless, non-relativistic bosons of mass m moving in a space of spatial 
■ coordinates of dimension d and stored either in a trapping potential U (r) or in a cubic box of 
' , ' size L with periodic boundary conditions. In this subsection the bosons are not interacting. The 
(— ( grand canonical Hamiltonian has then the following expression in second quantized formalism: 

o. /. 

O ■ H= dr^p\r){ho-Mr). (1.1) 

> ■ 

k>( _ jj, is the chemical potential, and the differential operator Hq includes the kinetic energy operator 
j_j ' and the trapping potential U (r) : 

ho = -^A + U{v) (1.2) 

where A is the Laplacian operator in dimension d. Note that [/ = in the spatially homoge- 
neous case of a cubic box. The field operator ^(r) annihilates a boson in point r and obeys 
the usual bosonic commutation relations: 



[V3(r),^(r')]-0, (1.3) 



(1.4) 
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It is convenient to expand the field operator on the orthonormal basis of the eigenmodes 
00, (r) of Hq with eigenenergy e„: 

V-W =^0<,(r)a„ (1.5) 



= eQ0Q(r) (1.6) 



where Oq, annihilates a particle in mode 0^. The grand canonical Hamiltonian is then a sum 
of decoupled harmonic oscillators: 

-ff = ^(ca - /^)alaQ. (1.7) 

a 

1.1.2 In thermal equilibrium 

We consider for convenience that the gas is in thermal equilibrium in the grand canonical 
ensemble so that the density operator of the system reads 

p^E-'e~^" (1.8) 

where /3 = {kBT)~^ is inversely proportional to the temperature T. As the canonical ensemble 
is more realistic, since no reservoir of particles is present in usual experiments, one has to check 
for each particular case at hand that the two ensembles are almost equivalent, a point that we 
will address in ll.1.41 

For the ideal Bose gas the density operator is a Gaussian in the field components so that 
Wick's theorem, an elementary proof of which is given in the appendix allows to express 
the expectation value 

(6) EE Tr[6/5] (1.9) 
of an arbitrary observable O in terms of the mean occupation numbers of the eigenmodes: 

Ua EE (alaa). (1.10) 

These occupation numbers are given by the Bose formula, also rederived in the appendix IXI 

^ (1.11) 



Assuming here for convenience that the ground state energy of /iq vanishes, which amounts to 
shifting the chemical potential, we see that the positivity of Ua imposes the following range of 
variation of /i: 

-oo<^<0. (1.12) 

In what follows one considers a temperature regime significantly different from the zero 
temperature case, that is one assumes that there is a large number of eigenmodes of energy 
less than ksT. In the spatially homogeneous case this implies that the size of the box is larger 
than the thermal de Broglie wavelength: 

L>A. (1.13) 

One recalls the various regimes then encountered for fj, increasing from — oo to zero, in the 
present case of a non-degenerate ground mode: 
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• ^ kBT: non-degenerate regime. All the occupation numbers are much smaller 
than unity and are well approximated by the Boltzmann formula: 

n„ ~ e-^(^"^^) « 1. (1.14) 

Summing this expression over a relates the chemical potential to the mean number N 
of particles. In the case of the cubic box, and if one replaces the sum by an integral, as 
allowed by condition H1.13|l . this leads to 

e^'^-pA'^ (1.15) 

where p — N/ L'^ is the mean density and A is the thermal de Broglie wavelength 

so that we recover the condition for a non-degenerate regime 

pA'^^l. (1.17) 

• degenerate regime: ^ fc^T so that the ground mode of the gas has a large occu- 
pation number: 

no ^ ^ » 1 1.18 

ImI 

where we have expanded the exponential in the Bose law to first order in its argument. In 
this case several modes have a large occupation number, as soon as <C fc^T, because 
condition H1.13|l holds. 

• the ground mode is more populated than any other mode: this regime is reached 
when the ground mode is more populated that the first excited mode: 

knT knT 

no ^ » ~ , , (1.19) 

which results in 

I^Kei- (1-20) 

It is the right time to recall the phenomenon of saturation of the total occupation of the 
excited modes, a key consequence of the Bose law: for a given temperature the maximal 
value of the excited modes population is 

^' ^ E < ^ E ^^^l (1-21) 

where one used the fact that each occupation number, being an increasing function of 
the chemical potential, is bounded from above by its value for p, — 0. An important 
consequence of (jl.2()|l is that the saturation of the excited modes population is quasi 
reached, as fi is negligible compared to all Sa in the sum defining A^4ax- 



(1.22) 
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• Bose-Einstein condensation: the ground mode has an occupation number on the 
order of the total population of the excited levels. One has then 

no ~ <ax- (1-23) 

One recovers in particular the condition to reach Bose-Einstein condensation applicable 
for a finite size system, that is in the absence of thermodynamical limit [1], 

N > Af^ax (1-24) 

and the typical value of the chemical potential when a clear condensate is formed, as- 
suming no ^ N [2]: 

ImI ~ (1-25) 



N 



1.1.3 3D vs ID in the thermodynamical limit 



We estimate the value of N^^^^ in a box of length much larger than the thermal de Broglie 
wavelength A: 

^^^^ = 5 exp(/3n^fcV2m)-l ^^'^^^ 

where the components of the wavevector k are integer multiple of 2-k/L. 

In 3D, the condition X allows to approximate the sum by an integral: 

N^,^c^am{j)' (1-27) 

where C is the Riemann Zeta function. One finds that the contribution to the integral is indeed 
dominated by wavevectors such that h^k^/2m ~ fc^T since the density of states cuts the low 
k contribution and the Bose law cuts the high k contribution. So in 3D in the thermodynamic 
limit, the gas is essential cither in the non-degenerate regime or in the Bose condensed regime. 

In ID one cannot replace the sum over k by an integral without getting an infrared divergent 
result for N^^: this is an indication that the low k dominate. The proper approximation to 
get the leading order contribution in L/X is to replace the Bose law by its low k behaviour: 

^-^-^^^^^ = ^1^ ^'-'^^ 

where we used Y^n=i = 7r^/6. So in ID in the thermodynamic limit, one can reach an 
arbitrarily high degeneracy pX without forming a Bose-Einstein condensate. 

1.1.4 Correlation functions of the ID Bose gas in the thermodynamic limit 
The first order correlation function of the field operator is defined in ID as 

gi(z) = (^t(^)^(O)). (1.29) 

It allows to determine the coherence length of the gas. In a ID box, gi is the Fourier transform 
of the momentum distribution: 
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where corresponds to the thermodynamic hmit. In the non-degenerate regime the momen- 
tum distribution reduces to the Boltzmann distribution so that gi is a Gaussian: 

In the degenerate regime, where <C fc^T, we replace the Bose law by its low energy approx- 
imation, as we did in the calculation of A^maxj which results in a Lorentzian approximation to 
the momentum distribution: 

rife ~ . (1.32) 

After integration over k one obtains the equation of state of the gas: 

and after Fourier transform, one obtains an exponential function for g\: 

5i(^)~pe-l^l/'= (1.34) 

with the coherence length 

lc = ^^X. (1.35) 

Note that in 3D, when a Bose-Einstein condensate is present, g\ tends at infinity to the 
condensate density po- 

The second order correlation function of the field to be considered here is 

g^{z) = (V:t(^)^t(o)v3(0)v3(z)). (1.36) 

Note that it is simply the spatial correlation function of the field intensity (that is the density) 
for 2; ^ 0. It is also proportional to the probability density, when measuring the positions of 
the particles, of detecting the first particle in and the second one in z, as is most apparent 
in first quantized form 

52(z) = ^(<5(f, - z)b{zi)) (1.37) 

where zi is the operator giving the position of particle i along z. 

For an ideal Bose gas at thermal equilibrium in the grand canonical ensemble, the second 
order correlation function is readily expressed in terms of g\ through Wick's theorem: 

g^{z) = ^ g\{z) (1.38) 

where we used the fact that g\ is real. This identity illustrates the phenomenon of bosonic 
bunching in real space: gi is maximum in 2; = where it reaches twice its asymptotic value 
. This implies strong density fluctuations correlated over a spatial scale on the order of the 
coherent length l^, a very negative fact from the atom laser perspective. 

What happens in 3D? The same identity for g-^, holds, as a consequence of Wick's theorem. 
When a condensate is present, one then finds also large density fluctuations at an arbitrarily 
large scale: 52 (+00) = + Pq> p^. Fortunately, this prediction is not correct physically and 
is an artifact of the anomalously large fluctuations of the number of condensate particles in 
the grand canonical ensemble for an ideal Bose gas. Wick's theorem implies 

{afal) = 2nl (1.39) 
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where ao annihilates one particle in the condensate mode, so that the standard deviation of 
the number of condensate particles is predicted to be equal to its mean value Uq, which is 
not correct physically! As a consequence the fluctuations of the total number of particles N 
becomes also anomalously large, A ~ ( A^) , and the equivalence between the grand canonical 
ensemble and the canonical ensemble is lost. This problem is well known for the calculation of 
the fluctuations of the Bose condensed ideal Bose gas [3]. An easy way to avoid this problem 
is to use the canonical ensemble, plus a number conserving Bogoliubov type approximation, 
see section 7.8 of [1] or [4]. 

In ID, when can we use the grand canonical ensemble to calculate fluctuations? In a box of 
size L, the variance of the total number of particles is related to the mean number of particles 
and to g2 by 

YarN^{N)+ dz / dz [g2{z - z') - p% (1.40) 
Jo Ja 

When no condensate is present in ID, L ^ Ic so that, using Wick's theorem and the expression 
(|1.34|l for gi, one gets 

which legitimates the use of the grand canonical ensemble. 



1.2 Construction of a classical field model 

The model of A^ bosons in ID interacting with a Dirac 6 potential is exactly solvable, see [5-7] 
for the case of repulsive interactions and [8] for the attractive case. However the calculation 
of the correlation functions 171 and f/2 is difficult, if one wants to access the full position 
dependence. A notable exception is the calculation of 32(0), which is linked to the free energy 
thanks to the Hellman-Feynman theorem and therefore requires only a knowledge of the energy 
spectrum [9]. We present here the classical field version of this model, which is also exactly 
solvable and in a much easier way, as done in [10, 11] . We restrict here to the repulsive case. The 
attractive case is indeed physically very different, as solitons can be formed and the canonical 
ensemble, rather than the grand canonical one, has to be used, which makes the classical field 
model more difhcult to handle. 

The intuitive idea to construct the classical field model is to replace the quantum field 
operator ^^{z) by a complex field tp{z). The Hamiltonian for the quantum field 



H ^ dz 



— (a,^t)(a^^) + ^^t2^2_^^t^ 

2m 2 



(1.42) 



is then replaced by the energy functional of the classical field 



Em] 



dz 



2m 



dz 



(1.43) 



As a consequence the thermal density operator for the quantum field is replaced by a probability 
distribution for the classical field: 



-1 

class ^ 



(1.44) 



The first order correlation function gi is then given by the average of ip*{z)i}'{0) over Pciass, 
and g2 is given by the expectation value of |'0|^(z)|-!/)|^(O). 
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A more systematic way of constructing this classical field model is to use the Glauber-P 
representation of the density operator, a common tool of quantum optics [12], defined by the 
identity: 

p = / V^P[{m{i>}){W\ (1-45) 



where the integral is a functional integral over all values of the complex field ^ and where we 
have introduced the Glauber coherent state of the matter field: 

IM) =e"/l'^l'/'e/''"'^^"^'^'^"^|0). (1.46) 

We recall that this coherent state, parameterized by the classical field 'ip{z), is an eigenstate of 
'ip{z) with the eigenvalue ^{z), for all z. We also note that the distribution P is not necessarily 
positive, nor even a regular function as it can be a distribution. The expectation value of 
normally ordered expressions, that is with all the on the left and all the V^'s on the right, 
are then exactly related to expectation values of classical field products with respect to the 
distribution P. For example 

j Vi^P[{^j}]r{zmO). (1.47) 

The classical field model is therefore defined by the substitution P Pciass- 

The effect of this substitution is well monitored in the case of the ideal Bose gas. The 
Glauber P distribution of the quantum field at thermal equilibrium is a Gaussian: 

P[{V'}] cx e"Sj"'=l'/"'= (1.4S) 

where the sum is taken over the eigenmodes of the one-body Hamiltonian, that are plane waves 
in the homogeneous case, ak is the amplitude of the complex field ip on the mode k and nk is 
the mean occupation number given by the Bose law [13]. The field distribution Pciass in the 
classical field model is also a Gaussian and is obtained by the substitution 

rik - nf- - (1.49) 
efc - A* 

where e^, here equal to h?k'^ /2m, is the eigenenergy of the mode k. One then recovers the 
equipartition theorem 

(efc - ^)n^.'-^^ = fcsT (1.50) 

as Boltzmann would have written it for a complex classical field at thermal equilibrium. We 
actually already used this expression for the occupation numbers, as a low energy approxima- 
tion to the Bose law in the degenerate regime, see 11.1.31 The classical field formula is indeed 
close to the Bose law for modes with a large occupation number: these modes of the field have 
almost negligible quantum fluctuations [14]. 

When is the classical field model expected to give predictions close to the full quantum 
theory? The answer depends on the observables and on the dimensionality of space. One has 
to remember first that the classical field model is subject to divergences in the absence of an 
energy cut-off, the most famous illustration being the black body catastrophe of nineteenth 
century physics. One should first introduce an energy cut-off emax on the energy of the modes, 
typically on the order of fc^T or larger. Some observables depend on this cut-off, like the mean 
energy of the gas, so they cannot be calculated with precision in the classical field model. Other 
observables, like the functions gi and §2 in ID (and this is a nice feature of ID), do not depend 
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asymptotically on the cut-off energy, as we have seen already for the ideal Bose gas and as is 
also the case in presence of interactions. In the ideal Bose gas case, the validity condition of the 
classical field approximation for the calculation of 51,(72 is, as we have seen, ksT ^ that 
is the gas has to be strongly degenerate. The same condition holds in the weakly interacting 
case, though with a different equation of state [11]: a more detailed discussion is postponed 
to the end of HI. 3. 51 Finally we note that the shift in the 3D critical temperature of a Bose 
gas due to interactions was calculated recently with a very high precision using a classical field 
model [15, 16]. 



1.3 Solution of the classical field model 

1.3.1 The problem to be solved 

One has to be able to calculate functional integrals of the type 

/I?V.V>*(z)7/;(0)e-^^M 



(1.51) 



where the complex field ip satisfies the periodic boundary conditions "0(0) — i'iL) [17]. The 
technique to implement this periodicity is to integrate separately over all the possible values 
V'o of V' in z = and over all possible closed paths starting in ■(/'o (for z = 0) and ending in 
•00 (for z — L). In the same spirit, in the numerator giving gi, one treats also separately the 
intermediate point of coordinate z, by integrating over the complex value ipz of the field in z. 
As a consequence 



r>vV'*(^)V'(0)e-'3^['^i 



i'(0)=ipo 



V(i)=V'o 



^(z)=i/'o 



^*(z)0(O) (1.52) 



where / df/'o stands for the integral over the real part and the imaginary part of 00 from —00 
to +00. One is then left with the calculation of functional integrals of the type 



(1.53) 



that is over possible complex fields defined over the spatial interval {zi,zf) with fixed values 
ipij at the end points of the interval. The energy of a field configuration over this interval is 
the spatial integral of the energy density 



Ez 



f * dz 




dip 


J Zi 


2m 


dz 



9 



(1.54) 



1.3.2 Reminder: Feynman's formulation of quantum mechanics 

Consider a quantum particle of mass M moving in 2D in a static potential V{x,y). The 
corresponding Hamiltonian is 



n 



Px 



2M 



Vix,y) 



(1.55) 
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where px,Py,x,y arc the momentum and the position operators along each direction of space. 
The so-called Feynman propagator in real time can be written in terms of a functional integral 
over two real valued paths x{t), y{t): 



(ar/,t//|e-«(*/-*')/'^|ari,2/,) 
where the Lagrangian is 



Vx 



x{ti)—Xi 



Vye 

y(ti)=yi 



£=-M 
2 



dt I 



dt 



V[x{t),y{t)]. 



(1.56) 



(1.57) 



We shall need this formulation of quantum mechanics for the propagator in imaginary time. 
Inside the integrals one performs the substitutions 



itf 
iti 
it 
idt 

where all the 'times' r's are real, to obtain 



n 

T 

dT 



fx(Tf)=Xf fV(Tf)=yf -T'dTSIn 

{xf,yf\e-^^^f-^'^/^\xi,yi)= Vx Vye->n 



rvi'rf)=Vf 



j;{Ti)—x 



which involves now the energy rather than the Lagrangian; 

2 / , \ 2' 



£ = -M 
2 



\dT J 



dT 



V[x{T),y{T)]. 



(1.58) 
(1.59) 
(1.60) 
(1.61) 



(1.62) 



(1.63) 



1.3.3 Quantum reformulation 

We map the ID classical field model at thermal equilibrium to a fictitious quantum propagator 
in imaginary time of a particle moving in a two-dimensional space, with the help of the following 
correspondence table. 



classical field model 


quantum analogy 


path integral 
abscissa z 

Im'(/'(2;) 


quantum propagator in imaginary time 
'time' T 
x{t) 

y{T) 



This correspondence is exact provided that j3 times the energy density of the classical field 
transforms into the energy £ over %, which imposes the mass M of the fictitious particle 

M = — — 1.64 
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and the fictitious external potential of the 2D quantum problem 



(1.65) 



The numerator of gi , given in H1.52|l , is transformed as follows with the quantum mechanical 
analogy: 



dtpo / dip. 



V'(z)=iAo 



Tr 



(1.66) 



The denominator of gi is then simply the trace of exp(— LTi/Ti). 

The calculation of the trace is performed in practice in the eigenbasis of TC. As this Hamilto- 
nian is rotationally invariant, one classifies its eigenstates with the angular momentum quantum 
number integer from — cx3 to +00, and the radial quantum number n, integer from to +00. 
The corresponding normalized eigenvector is denoted {(pl,). For a given I, the eigenenergies are 
sorted in increasing order starting from n — 0. The absolute ground state is obtained for Z = 
since a non-vanishing angular momentum gives rise to a centrifugal energy term that can only 
increase the energy. A great simplification occurs in the thermodynamical limit, that is when 



i(ecxc - e'„-°)/h » 1 



(1.67) 



where Ccxc is the lowest energy eigenvalue above the ground state energy ep"^". In this case 
one can approximate the imaginary time evolution operator by restricting to the lowest energy 
mode(s), e.g. 



-LH/h 



In the thermodynamic limit, one therefore obtains the expressions 



9i{z) 



bf\x-^y)e-^^-^o^°y\x + ^v)\cpf^) 



g,{z) = (0ri(x^+y^)e 



(i'+r)i0o 



i=0\ 



(1.68) 

(1.69) 
(1.70) 



1.3.4 Results 

We briefly present the predictions of the classical field model, more details being given in [11]. 

A first point is to realize that, with the proper system of units, the state of the classical field 
in the thermodynamic limit is controlled by a single dimensionless parameter. If one expresses 
tp in units of the square root of the mean density p and the length z in units of the coherence 
length of the ideal Bose gas, 

z^4^Z^(^Z (1.71) 

where Z is dimensionless, and if one realizes that the chemical potential /i is fixed by the 
condition 

{\^?)=P (1-72) 



WILL BE SET BY THE PUBLISHER 



Prl-11 



one can check that the classical field energy functional E in H1.43|l depends on the single 
parameter 



X ■ 



P9 



(1.73) 



Our concern for atom laser applications was the large intensity fluctuation of the ideal Bose 
gas. In the classical field model in the thermodynamic limit, the function 172 is a decreasing 



function of z € (0, +00), since TL 
we introduce the contrast 



is positive. It reaches its maximum value in z = so 



C 



(1.74) 



This contrast is plotted as function of x in figure ^ For x = 0, one recovers the ideal Bose 
gas value C = 2. When x increases, there is a very sharp decrease of x, followed by a power 
low tail, C tending slowly to unity when x tends to +00. The curve shows that repulsive 
interactions efficiently reduce the density fluctuations of the gas. 




Fig. 1. In the classical field model, contrast C of the local density fiuctuations of a ID Bose gas as 
function of the temperature dependent interaction strength x defined in Ea. H1.73|l . 



What is the physical origin of the very fast decrease of C in the vicinity of x = 0? This 
sharp crossover from strong to weak density fluctuations corresponds to a change of shape of 
the fictitious potential V . When x < Xc — 0.28, the chemical potential is negative so that V 
is minimum in x = y = see figure and the field remains trapped around = 0, with 
large intensity fluctuations. When x > Xc, M is positive so that y is a Mexican hat potential, 
see figure Eb, and the field is trapped around a non- vanishing value of its modulus, with weak 
density fluctuations. 

The large z behaviour of the correlation functions is also easy to access. In H1.69|l . x + iy 
maps \4>(r'^) onto a Z = 1 state so it is sufficient to inject a closure relation on the eigenstates 
of Ti within the I = 1 manifold. When z +00 one can keep the contribution of the mode 
n = 0,1 — 1 so that gi is proved to vanish exponentially with z: 



j=i 



j=o 



9i[ 



aie 



with 



(1.75) 
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(a) 





, 2 2,1/2 ,1 x1/2i 

(x+y) [(p) ] 



, 2 2,1/2 r, ,1/2, 

(x+y) [(p) ] 



Fig. 2. In the quantum mechanics equivalent to the classical field problem, potential V{x,y) seen by 
the quantum mechanical particle, (a) For a low value of x = 0.1 and (b) for a large value of x = 20. 
The units on the axes are such that V depends only on the parameter x defined in Ea. (ll.73|l . 



where ai = \ {(p\f^\x + iy\(j)'ij^'^)\'^ does not depend on z. We may therefore define the coherence 
length of the gas as For 52 the closure relation is in the I = manifold and one keeps 

the contribution of the ground and the first excited mode in the large z limit: 

g^iz) ^ gliO) + a2e---' with ^2 = ~ (1.76) 

where 02 = + y^\4'(r'^)\'^ is a constant. The correlation length of the gas may be 

defined as 1 /k2 ■ 

The values of ki, K2 are plotted in figureOl For g = one recovers the ideal Bose gas results 
of ll.1.41 For repulsive interactions, the coherence length is slightly increased, by a factor of up 
to 2. The correlation length strongly decreases and its ratio with the coherence length tends 
to zero in the large x hmit: the density fluctuations not only decrease but also take place on 
a length scale much smaller than the coherence length of the atom laser. 

Analytical results can be obtained for ki and K2 in the large x hmit by using the quantum 
mechanical analogy. Writing the wavefunction with angular momentum I as 

= (1.77) 



where r and 9 are polar coordinates in the x, y plane, one has to solve Schrodinger's equation 

-^^I"{r) + U{r)f{r)=ef{r) (1.78) 



with the effective potential 



U{r)^V{r) + -—^{l'-\/A) (1.79) 



and the boundary condition /(O) = 0. Note the presence of a centripetal potential for I = 0, 
which does not support any bound state, a known peculiarity of the 2D quantum motion. In 
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Fig. 3. In the one-dimensional classical field model, the inverse field coherence length ki and the 
inverse intensity correlation length K2 as function of the temperature dependent interaction parameter 
X defined in Eg. lirT^ . 



the large x limit one finds the minimum of U{r) in ro > and one approximates U{r) around 
To by a series expansion in powers of r — tq. To leading order, U(r) is then approximated 
by a quadratic potential centered in rg, then cubic corrections can be taken into account 
perturbatively, and so on. One has to be careful in this type of expansion to collect all the 
terms of a given order. For example, the calculation of the mean density proceeds through the 
identity 

(r^) =rl + 2{r - ro)ro + ((r - ro)^) (1.80) 

where the expectation value is taken in One has then to take care of the fact that 

the last two terms are of the same leading order, the last term being non-zero already for the 
harmonic approximation to U{r) whereas the other term (r — rp) is first non-zero when the 
cubic distortion of U{r) is included. One obtains in the large x limit: 



C 



P9 

TT 

1 

2 
1 



1 



2^1/2 



deal 



1 



1/2 



X 



1/2 



(1.81) 
(1.82) 

(1.83) 
(1.84) 



where ^ is the so-called healing length. The equation of state H1.81|l exactly coincides with 
the Gross-Pitaevskii equation in the large x limit, and ^ also naturally arises from the Gross- 
Pitaevskii equation. This may be surprising at first sight, since the Gross-Pitaevskii equation 
holds for the condensate wavefunction, whereas there is no true condensate in ID in the ther- 
modynamic limit. This fact will be recover and hopefully will become clear in section |21 of 
this lecture. We note that x is proportional to {Ic/S,)^, so that the regime of weak density 
fluctuations corresponds to Ic 3> ^. 
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1.3.5 A Bogoliubov approach for x ^ 1 



A limitation of the previous analytical expansion in the large x liniit is that it is difficult to 
transpose to the quantum problem. Is there another approach available that would apply both 
to the classical field and the quantum field? 

One could think of applying the Bogoliubov approach to the classical field model. A problem 
appears in the large L limit since no condensate is present, as we now show. In the classical 
field version of the Bogoliubov approach, one splits the field as 



1 



Ll/2 



(1.85) 



where ip^ is the component of the field on the plane waves with non-zero momentum and is 
supposed to be much smaller than the component on the A; = wave, oc y^No/L, so that one 
can neglect all the terms of degree > 2 in tp± in the energy function (|1.43(l . The resulting 
quadratic function of the field can then be diagonalized by the Bogoliubov transformation. As 
shown in the lecture by Gora Shlyapnikov in this school, one then has 



dz{\i.^{z)\')=Y,{Ul + V,') 



bog 



(1- 



where Uk, Vk are the dimensionless amplitudes of the Bogoliubov modes on the plane waves, 
normalized such that — 1, and e^°^ is the usual Bogoliubov spectrum. This sum is 

dominated by the low momenta: 



m 



As a consequence the relative mean weight of the excited plane waves is 



I Hi 



mkBT 1 



KlL 

~6~ 



(1.87) 



(1- 



where the inverse coherence length is given by (|1.82|l . We see that the basic assumption of 
Bogoliubov theory fails when L ^ kJ"^, that is in the absence of a condensate. 

The appropriate small parameter in the large x domain is not the non-condensed fraction 
but the relative density fluctuations. This gives the idea of using the intensity-phase represen- 
tation of the field, an idea extensively developed in [18]. One writes 



one recalculates the energy functional of (|1.43|l : 



Jo 



2m 



dz 



VP 



9 2 

2^ 



MP 



One then introduces the density deviation 5p such that 

p{z) = po + 5p{z) 



(1- 



(1.90) 



(1.91) 



1/2 

and po is the constant such that the uniform field V'o — Pq has the absolute minimum of 
energy: 

M = 5Po. (1.92) 
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Finally one quadratizes the energy functional in terms of 5p and 9. Remarkably this does not 
require weak phase fluctuations, only low density fluctuations, so that 



(1.93) 



and this is a crucial point in the absence of a condensate! This leads to the quadratic energy 
functional 



ct 



dz 



(1.94) 



This energy functional may be put in normal form by a Bogoliubov transformation. If one 
imposes periodic boundary conditions for the phase 9, one obtains exactly the usual Bogoliubov 
spectrum. This allows to recover the large x analytical results of ll.3.41 [19]. 

This modified Bogoliubov approach can tell when the classical field model gives predictions 
close to the full quantum theory for % 3> 1: 

• the shortest length scale obtained from the gi , g2 functions is the correlation length 
~ ^, corresponding to a Bogoliubov energy of the order of p. The corresponding 
Bogoliubov modes should have a large occupation number for the quantum field to be 
approximated by a classical field, which imposes 



ksT :s> p ~ pg. 



(1.95) 



• the gas should be degenerate. Since (pA)^/27r ~ x^sT/pg, the condition pA 1 is 
automatically satisfied for a large x when 11.95|l holds. 

• Also a classical field model is realistic in the weakly interacting regime only: strong 
correlations between the particles are absent in a classical field model . . . since there is 
no such a thing like a particle! This is apparent from the equation of state p = pg, which 
makes sense quantum mechanically in the weakly interacting regime only, that is when 

3> 1. Since pS, = X^^^f^sT/ pg, pS, is automatically big in the large x regime when 
Cnm|l holds. 



2 EXTENSION OF BOGOLIUBOV METHOD TO QUASI-CONDENSATES IN THE 
WEAKLY INTERACTING REGIME 

The ID classical field model of the previous section has identified an interesting regime for 
guided atom optics and atom lasers, a degenerate regime where the density fluctuations are 
strongly reduced by the repulsive interactions between the particles and where the coherence 
length is much larger than the thermal de Broglie wavelength. Such a state of the gas exists even 
if there is no Bose-Einstein condensate; it is called a 'quasi-condensate', and was extensively 
studied in [18] for the uniform case (see the lecture of Gora Shlyapnikov for the trapped case). 

The goal of this section is to present a very simple but accurate theoretical frame, simpler 
than the one of [18], to study the quasicondensates in the degenerate and weakly interacting 
regime. Degenerate means that 

pX'^ > 1 (2.1) 

where p is the gas density, A is the thermal de Broglie wavelength and d is the dimension of 
space. Weakly interacting means 

» 1 (2.2) 
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where the healing length ^ is related to the chemical potential of the gas by 



M = 



(2.3) 



Our theoretical frame is simply the Bogoliubov approach but in the density-phase representa- 
tion of the quantum field. In case where a condensate is present, it gives results equivalent to 
the usual Bogoliubov approach. 



2.1 Construction of an appropriate model 

A difficulty often appearing in theoretical treatments of the quasi-condensate problem is the 
appearance of divergences, an infrared divergency in ID, an infrared and an ultraviolet diver- 
gency in 2D. We construct here with some care a model Hamiltonian that, when combined with 
a systematic expansion in powers of small parameters, will allow to avoid all these divergencies. 



2.1.1 General idea of the theory of quasi-condensates 

This general idea is for example inspired by the Bogoliubov approach developed in the previous 
section for the classical field model. One identifies as the small parameter the relative density 
fluctuations of the gas and one performs a systematic expansion in this small parameter. So 
one flrst writes the field operator as 



V'(r)=e'^«V^ (2.4) 

where the operator giving the density is 

p(r)=Vi+(r)Vi(r) (2.5) 

and where 9 is the mythical phase operator in point r. Theses two operators are expected to 
be canonically conjugated: 

[p{r),e{r')]=i6{r-r'). (2.6) 
One then splits the operator giving the density as 



p{r)=po+Sp{r) (2.7) 
where po is a number, and one quadratizes the Hamiltonian 



H = J ^t/jQ^ + |^t2^2 (2.8) 

in dp under the condition of weak density fluctuations, that is of a small variance of p: 

Varp = _ (p)2 « p2 (2.9) 



where p is the mean density. 
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2.1.2 Problems with the general idea 

Whereas the general idea works fine with a classical field model, it is plagued by three major 
problems for the quantum field problem. 

Problem 1: the representation of the atomic interaction potential by g times a Dirac 
distribution is perfectly fine in ID, but leads to a mathematically ill-defined quantum problem 
in 2D and in 3D and gives rise to ultraviolet divergencies. 

Problem 2: The variance of p is infinite, so that condition 1)2.9(1 is not satisfied. The 
second moment of p can indeed be expressed in terms of the g2 function by using the bosonic 
commutation relations for -0: 

=5(0)^ + 32(0) =+(X3. (2.10) 

Note that the function g2 is bounded from above in any realistic model. This divergence comes 
from 'quantum fluctuations' of the field, that is from the fact that is not normally ordered 
with respect to -0. 

We conclude that condition (|2.9() is not the proper physical way of defining the regime 
of weak density fluctuations. Everything becomes clear if one considers the statistics of the 
number of particles in a finite volume of space, e.g. a box of length / and volume The 
operator giving the number of particles in the box B is 

<i''r7/'t(r)v3(r). (2.11) 

B 

This operator has now finite relative mean squared fluctuations: 



Yarn 1 f d'^r f d'^r' 



Id I Id 



52(r-r') 



(2.12) 



The first term on the right hand side of this expression is the quantum term leading to the 
infinite variance of p: it diverges in the limit I —>■ 0. In the treatment to follow, we will consider 
a size of the box large enough so that 

(n) = pl'^ > 1. (2.13) 

In the regime of weak density fluctuations the second term in the right hand side of H2.12|l is 
also small. If one can find a length I smaller than the spatial scale of variation of g2, but still 
satisfying pl'^ ^ 1, the regime of weak density fiuctuations corresponds to 

|^-1|«1. (2.14) 
P 

Problem 3: there exists no hermitian phase operator satisfying (|2.6|) . We produce a 
simple proof ad ahsurdum [22]. Imagine that there exists a Hermitian operator satisfying the 
commutation relation l|2.6|l . Then the following operator T(a) is unitary: 

f{a) = e^*"^^'"") (2.15) 

where a is real number and rp is a fixed point in space. Let us calculate the corresponding 
unitary transform of n: 

n{a) = f\a)nf{a). (2.16) 
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Fig. 4. Two possible sliapes for tlie probability distribution ttk of the number n of particles in a box 
B of volume l'^. (a) For a quasi-condensate with an appropriate choice of I: the most probable value 
is > 0, 7r„ is approximately Gaussian with a standard deviation smaller than the mean value, (b) For 
an ideal Bose gas in the absence of Bose condensate: the most probable value is n = 0, 7r„ is essential 
an exponential function of n. The mean value of n is 8 in both examples. 

Assume that the point Tq is inside the box B. Then the commutation relation (|2.6|l implies 

[^(ro),n] = (2.17) 
As a consequence, n(a) obeys the differential equation 

■^n{a) = e'"^('"«'z[0"(ro), nle-*"^^'"") = 1 (2.18) 
da 

with the 'initial' condition h{0) = n. This finally leads to 

h{a)=h + a. (2.19) 

We then have the following disaster: if lip) is an eigenvector of n with eigenvalue n, T{a)\^) is 
an eigenvector of n with eigenvalue n + a. As a can be any real number, this contradicts two 
well established properties of the spectrum of n, the discreteness and the positivity. 

Violation of discreteness is somehow attenuated by the fact that in practical physical calcu- 
lations, only integer values of a appear: in the Hamiltonian and in the gi function for example, 
one has only a = 0, +1 or —1. Violation of positivity is in general a real problem. It becomes 
in practice a minor problem if the physical state of the system is such that 

(n) > 1 and Varn<(n)2 (2.20) 

and such that the probability of having or 1 particle in the box B is absolutely negligible. 
In other words the use of the commutation relation (|2.6(l is an acceptable approximation when 
the probability distribution 7r„ of the number of particles in the box B has a shape like in 
figure which is typical of a quasi-condensate, but not when 7r„ has a shape like in figure 
which is typical of the non-condensed ideal Bose gas. 

2.1.3 A solution to these three problems 

The introduction of a finite size box is a key ingredient of the previous subsection 12. 1.21 We 
therefore discretize the real space on a grid with a step a square grid in 2D and a cubic 
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grid in 3D. The position r now refers to the coordinates of the center of each cell. The field 
operator V'(r) annihilates a particle in the cell of center r and is normalized in such a way that 

[^(r),^t(r')] = ^ (2.21) 

to recover the usual continuous quantum field theory for I — s- 0. We use periodic boundary 
conditions with lengths that are integer multiple of I. 

Solution to problem 1: the interaction potential is modeled by a discrete delta potential 

nr-r') = |<5r,r' (2.22) 

where the choice of the so-called bare coupling constant go is discussed later. Note that, in 2D 
and 3D, go strongly depends on I in the limit Z — s- 0, which is a signature of the pathology of a 
Dirac delta potential in a continuous 2D or 3D space. 

The grand canonical model Hamiltonian in second quantized form is then 



(2.23) 



where we have introduced the discrete Laplacian A. In the initial stage of the Bogoliubov 
approach to come, it is convenient to represent this discrete Laplacian by 

A/(r). E fir + leMfiyle,)^2fir)^ ^^^^^^ 

j^d directions of space 

In this case our Hamiltonian becomes fully equivalent to the Bose Hubbard model with an 
on-site interaction U = go/l"^ and a tunneling amplitude t = —h^/2mP. The price to pay is 
that the dispersion relation of the kinetic energy term is 

^2 d 



^k" = ;^E(l-cos^^O (2.25) 



which coincides with the correct parabola only in the range kl ^ 1. In the final use of the results 
derived by the Bogoliubov approach, in particular in a numerical treatment for the spatially 
inhomogeneous case, it is therefore more accurate to define A in Fourier space, simply by the 
requirement that the plane wave of wavector k is an eigenstate with eigenvalue H^k^ /2m, each 
component kj of the wavevector being restricted to the first Brillouin zone [—n/l, 

Solution to problem 2: the variance of the density is finite. In the discrete model, the 
operator giving the density is p{r) = ip'^{r)ip{r) and has the physical meaning of being equal to 
l~'^n where h is the operator giving the number of particles in the considered cell. The same 
phenomenon as in i j2. 1.21 takes place, the variance of p is finite: 

Varp=^ + [(^t2^')-p'] (2.26) 

where p is the mean density in the considered cell. We assume that both terms in the right 
hand side of (|2.2t)|l are much smaller than p^: 

pl'^ > 1 (2.27) 

(2.28) 
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The first condition ensures that the cell has a large occupation, and the second one ensures 
that the density fluctuations are weak. 

Solution to problem 3: approximate construction of the phase operator. Following 
[23,24] one can introduce the exact writing 

^ = A^fp (2.29) 

where A decreases the number of particles in the cell by one with an amplitude one rather 
than ^Jn: 

A\n) = \n-l) if n > (2.30) 

= otherwise. (2-31) 

The following identities are exact: 

lit = 1 (2.32) 

= 1-|0)(0| (2.33) 

where |0) is the vacuum state of the considered cell. This reveals that A is not unitary. However 
if the probability of having no particle in the cell is truly negligible in the state on which A acts, 
the overlap of this state with |0)(0| can be neglected and one can assume that A is unitary: 

iU ~ iit = 1. (2.34) 

One then can set 

A = e'^ (2.35) 

where 9 is Hermitian. The commutation relation (|2.6|l is then an acceptable approximation 
that ensures e.g. that 

itpi-/5-^. (2.36) 

We note that there exists rigorous definitions of the phase operator [25], but the resulting 
operator does not satisfy (|2.6() . which makes explicit calculations more difficult. 



2.1.4 How to choose the grid spacing I ? 

A first condition is that I is large enough so that pl''^ ^ 1. Note that the model Hamiltonian 
H2.23|l can be used outside this regime, this condition being useful only for the Bogoliubov 
approximation to come. 

A second condition is that I is small enough so that the energy cut-off '-^ fi^ /mP introduced 
by the grid does not change the physics. This obviously requires that the energy cut-off is 
larger than /cbT and than the chemical potential p. Equivalently 

l<\£. (2.37) 

where A is the thermal de Broglie wavelength and ^ is the healing length defined in (|2.3|) . 

These three conditions are compatible only in the degenerate ()2.1f) and weakly interacting 
H2.2|l regime to which our Bogoliubov approach is therefore restricted. 
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2.1.5 How to choose the couphng constant go ? 

The idea is the following. Consider the scattering of a plane wave on the discrete S potential on 
the grid, for a finite grid spacing I but of course in the case of an infinite quantization volume. 
Calculate the corresponding scattering amplitude and compare it to the exact amplitude for 
the true interaction potential in continuous space. Adjust go to have the same scattering 
amplitude in the low energy domain. 

Let us calculate the T matrix for two interacting particles in the discrete model. We take 
the version of the discrete Laplacian giving the correct parabolic spectrum, as discussed after 
H2.25|l . The center of mass motion can be separated from the relative motion so we consider 
the reduced Hamiltonian for the relative motion of the two particles on the grid: 

H= — + V (2.38) 

m 

with V — go5r,o/l'^- The eigenstates of the kinetic energy are plane waves with a wavevector 
k, with a wavefunction 

(r|k) = e^''"'. (2.39) 

As each component Tj is an integer multiple of Z, the component kj has a meaning modulo 
2-k/I and can be restricted to the first Brillouin zone [— tt/Z, 7r/Z[. We normalize the localized 
state vector |r) so as to recover the continuous theory for I — > 0: 

(r|r') ^ (2.40) 

In this case the potential V can be written as 

^ = go|r = 0)(r = 0| (2.41) 
and the following closure relation holds: 

(0i|k>(kHl (2.42) 

where D = [— tt//, 7r//['' is the first Brillouin zone of the reciprocal lattice. 
The T matrix at a given energy E is defined by [26]: 

T{E + iri) = V + VG{E + iri)V (2.43) 

where 77 ^ 0+ at the end of the calculation. G is the resolvent of the full Hamiltonian: 

G(z) = -\ (2.44) 

defined for any complex number z not belonging to the spectrum of H. The particular form of 
V leads to the simple expression for the matrix element of T in between two arbitrary plane 
waves [27]: 

(k\TgM{E + n])\k') = go + .gg(r = 0|G(S + M;)|r = 0) (2.45) 

so that only the matrix element of the resolvent in the state localized in the cell r = of the 
lattice matters! This matrix element is immediately deduced from the recursion relation 



G{z) = Goiz) + Go{z)VG{z) 



(2.46) 
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where Go is the resolvent of the kinetic energy operator. We finally obtain 
where 

m = {v = 0\G,{E + ir^)\v = Q)= lim / f'^,^^. ^2,,, • (2-48) 

Note that, due to the discrete delta nature of the interaction potential, the matrix element 
of Tgrid depends only on the energy, not on the wavevectors: the operator Tgrid is actually 
proportional to |r = 0)(r = 0|. 

To calculate 1{E) one may use the following identity for distributions: 

lim --^ = VV\- - i-jTS(X) (2.49) 
^^0+ X + ir] X ^ ' ^ ' 

where VV is the principal part. We restrict to positive energies E and we set E — fi'q^ jm. We 
are interested in the low energy limit E <C /mP that is g <C tt/L Why ? In the degenerate 
and weakly interacting regime the maximal value of g is 1/^ or 1/ A, which is much smaller 
than l/l. 

The imaginary part of I{E) is easy to calculate when q < n/l, in which case the support 
of the S distribution is inside the Brillouin zone D and Im I{E) is proportional to the density 

of states: 

TT) 

Im 1(E) = 5-g in 3D 

TfX 

—2 in 2D for g < tt/Z (2.50) 



An 

m 1 
2n^ q 



in ID 



The calculation of the real part of I{E) is more involved: 

ReI{E) = ^ [ -^VV^^ (2.51) 

so we shall restrict to the low energy limit. In 3D, the real part has a finite limit for q ^ 0: 
replacing q by zero, one obtains 



Re7(£; = 0) 



Tfi 

0.194... (2.53) 



The situation is dramatically different in 2D, where one gets a divergent expression if one 
replaces q by zero. The trick is then to split the integration domain D in the disk of radius 
tt/Z, over which the integral can be calculated exactly, and in the complementary domain, 
where q^ is a small perturbation of k'^ [28] : 

= ^ [ln(gZ/7r) -C + 0{{qlf)\ (2.54) 
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with the constant 

2 

C = In 2 - -G = 0.110025 . . . (2.55) 

TT 

involving Catalan's constant G — 0.915965594 .... In ID the exact calculation can be done: 



777 

Re/(£;) = ^In 



TT + ql 



TT — ql 



ml 

^,1^0 (2.56) 
TT^n 



The last step is to compare with the scattering amplitude of the true interaction potential 
in the continuous case. In 3D a low energy approximation can be used for the true T matrix: 

(k|rt,uc(g + »7?)|kO , '. , (2.57) 

a +^9 

where a is the s-wave scattering length, provided that kr^, k've, qr^ <^ 1, where r^. is the 
effective range of the true interaction potential. The maximal value of q, k, k' that one may 
expect to appear in the thermal state of the gas in the degenerate regime is the inverse of the 
mean interparticle separation p^^^, whatever the strength of the interactions: one has to check 
that prl <C 1, which is the case in present experiments even close to a Feshbach resonance [29]. 
In the weakly interacting regime, the situation is even more favorable as the typical q, k, k' is 
at most 1/^ or 1/A, much smaller than p^^"^. 

Using the expressions (|2.53|l and H2.50|l for I{E) and identifying Tgrid with Ttiuc leads to 
the identity 

1 m m f (Pk 1 



go ATTh^a Jd (27r)3 k^ 



(2.58) 



that is 



g, = p!L^ With i^3- 2.442 749... (2.59) 
1 — K^a/l 

Note the exact cancellation of the imaginary parts of the denominators of Tg^id and Ttruo- In 
the considered weakly interacting regime, condition (|2.2() implies \/~pa^ <C 1; as pl^ ^ 1 one 
finds that a <C / so that go is close to the usual coupling constant g = ATTh^a/m. Close, but 
not identical, and this plays an important role in suppressing potential ultraviolet divergences. 
In 2D, the low energy approximation for the 'true' T matrix is 

(k|Tt™e(i? + iriW) ■ , ^If^"^ . ,^ (2.60) 
ln(a2g/2) + 7 - itt/z 

where a2 is the 2D scattering length and 7 = 0.57721 ... is Euler's constant, see the lecture of 
Cora Shlyapnikov in this volume. In identifying with the low energy expression of Tg^id, one 
finds again that the imaginary parts of the denominators of Tgrid and Ttruo exactly cancel and 
one is left with: 

2TTh^ 1 



with a numerical constant 



where we used the amusing fact that —7 + 2G/tt = 5.91 ... 10 <C 1 



50 = 1 . (2.61) 

m \n{K2l/a2) 



K2 - ie-^+2GA ^ i (2.62) 

TT TT 

3 
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In ID one may model the continuous interaction potential by giS(zi — Z2), under validity 
conditions discussed in [30,31], in which case the T matrix is 

(fc|Ttruc(i? + ^ (2.63) 



where the ID scattering length is such that 



51 = — • (2.64) 
max 

The identification of T^^d and Tirue leads to 

= ~mi7 2T- (2.65) 

In the weakly interacting regime, go is very close to 171: / ^ ^, so that 1/ ai <C £,/ai oc 1/ ^fpal cx 
<C 1. Contrarily to the 3D case, the small diff'erence between go and gi does not play a 
significant role. 

2.2 Perturbative expansion of the model Hamiltonian 

In a regime where each cell has a negligible probability of being empty, we may introduce 
the phase operator having approximately the commutation relation (|2.6|l with the density and 
write the Hamiltonian as 

H = Hkin + i?pot (2.66) 



Hpot = 5]/'/5(r) 



C/(r)-A^+f (p(r)-^ 



(2.67) 



where we have introduced the notation 

/+, - /(r + /e,) (2.69) 

= /(r-/e,) (2.70) 

for each direction of space. 
2.2.1 Two small parameters 

The first small parameter expresses the weakness of the density fluctuations: one splits 

p(r) =po(r)+<5/3(r) (2.71) 

where po, the zeroth order approximation to the density, corresponds to a pure quasi-condensate. 
The small parameter is then 

.1 = M « 1. (2.72) 

The second small parameter expresses the smallness of the phase variation in between 
neighboring sites of the grid: 

£2 = l^grad^l < 1 (2.73) 
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where grad is the discrete gradient on the lattice. 

It can be checked at the end of aU calculations that by an appropriate choice of I one can 
achieve 

^ (2.74) 



ei ~ £2 



One may note that this value corresponds to the 'quantum' fluctuation part, that is the first 
term in the right-hand side of H2.26|l . An important consequence is that the two small param- 
eters may be considered as a single small parameter, in the expansion of the Haniiltonian. 



2.2.2 Quadratisation of the Hamiltonian 

Expansion of H in powers of ei and e2 is performed from the series expansions 

5p Sp'^ 6p'^ 



1 



2pI^' 



3/2 



ic 5/2 



(2.75) 



(2.76) 



To zeroth order in the small parameters ei and £2, density fluctuations are neglected, p being 
approximated by po, and the spatial variation of the phase operator is also neglected. This 
leads to the zeroth order approximation to the Hamiltonian, which is the following c-number: 



'pT 



1 



^7^A + C/(r)-f -goPo~M 
2m 2 



1/2 
Po • 



(2.77) 



In 3D, this is the equivalent on the lattice of the Gross-Pitaevskii energy functional, with the 
difference that the bare coupling constant t/o appears, rather than g. pQ is then naturally 
obtained by minimization of Hq, so that it obeys: 



' — A + U{r)+gapo-fi 
Zm 



1/2 n 

Po =0. 



(2.78) 



This equation naturally defines po as a function of fi. However it will reveal more convenient to 
parameterize the theory in terms of A'o, i.e. the total number of particles stored in the density 
profile Po'- 

No = Y,l''Po{r)- (2.79) 



We will therefore consider p, and po as functions of Nq: 

P = Pq{Nq) 



Po(r) 



Poir;No). 



(2.80) 
(2.81) 



As a consequence of Ea. (|2.78|l . one finds that the first order approximation Hi to the 
Hamiltonian exactly vanishes. The first non-trivial contribution is therefore a quadratic one: 



H2 = E2[po]+Y.l' 



5p 



5p 



SmpQ^'^ 



90 



2m 



^ Jpo(i-)po(i- + ?ej 



{e{v + le;)-e{v)f 



(2.82) 

(2.83) 
(2.84) 
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where the c-number energy functional E2[pa] is given by 



^2N = -f E 



Po 



Amp 



EE 

r j=l _ 



Po 



1/2 



Po 

P0.+] 



1/2 



(2.85) 



Remarkably, one finds that the Hamiltonian H2 is equivalent to the Hamiltonian of the U{1)- 
symmetry breaking Bogoliubov approach for the usual case of true condensates: one introduces 
the field 



m-^^+^pTimr). 

One can then check that it has bosonic commutation relations 



[B{v),B\v')] 



Sr,r> 
Id ■ 



(2.; 



(2.87) 



After some algebra, using the fact that Pq^^ solves the Gross-Pitaevskii type equation H2.78|l . 
one obtains the following exact rewriting of H2: 



h 



H2^l''Y.B^ [- — A + U + goPo- f^) B + goPo 



13^ B + -{B^ + B^^) 



(2.88) 



which has exactly the structure of the usual Bogoliubov Hamiltonian. Remarkably, the c- 
number energy E2 is exactly compensated by the contribution of commutators of 6p with 9. 
One then reuses the standard diagonalization of the Bogoliubov Hamiltonian: 



B(r) = -zQpl/\r) + PdNoPrir) +J2bsUsir) + blv:ir) 



{2.. 



as explained in [24,32,33]. Here the sum s is taken over the regular Bogoliubov eigenmodes 



( -hi^ + U - f^ + '^9opo gopo 

\ -goPo - (-|^A + C/-Ai + 2goPoj 

normalizable as 

E^'[l«.(r)P-K(r)P] =1. 



(2.90) 



(2.91) 



Thermodynamic stability, or equivalently the fact that po is a minimum of Hq, imposes that 
the corresponding eigenenergies £5 are positive [1]. The operators bs obey the usual bosonic 
commutation relations [Ssj^g/] — Sg.s'- The operators Q and P commute with and are 

conjugate quantum variables, 

[P,Q]^~i. (2.92) 

Q is a collective coordinate representing the quantum phase of the field, and, as shown in [28] , 
P gives the fluctuations of the total number of particles away from the total number of particles 
No contained in the density profile po: 



P ^ 1'Y.^pir) - p„ir)] ^ N - No. 



(2.93) 
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We note that our treatment is in the grand canonical ensemble but does not break the U{1)- 
symmetry of the Hamiltonian: the emergence of the operators P and Q is therefore a conse- 
quence of a non-fixed value of the total number of particles, rather than of symmetry breaking. 

Replacing B by the Bogoliubov modal decomposition finally leads to a modal expansion 
for the operators giving the density and the phase: 

^(r) = -Q + J20sir)k + e:ir)bl (2.94) 

S 

5p{v) = PdN,p^{v)+Y.5ps{v)hs+5p:{v)hl (2.95) 

s 

where 

' ^ O- 1/2/ N 

5ps (r) = (r) [us (r) -t- (r)] 

It also gives the normal form of the Hamiltonian H2, from which thermal averages can be 
evaluated easily: 

2 'dNo ^ ^ ^ ' 2dNo 
where the function po is defined in Eq. (|2.8U|) [34] . 



(2.96) 
(2.97) 



2.2.3 Why cubic terms in the Hamiltonian are required 

One could believe that the knowledge of H2 is sufficiently to calculate to a given order 0(e^ 2) 
the correction for any observable to the pure quasi-condensate approximation. However this 
is not true, the mean density being an obvious counter-example: in a thermal state density 
operator with the Hamiltonian H2, the mean value of (5/5 vanishes so that, for a given chemical 
potential, there is no correction to the mean density as compared to the pure quasi-condensate 
assumption. 

A similar phenomenon takes place in the quantum mechanical problem of a single particle 
strongly confined in a non-harmonic trapping potential V{x) in ID: a series expansion of the 
potential is performed around its minimum, 

Vix) ^Vo + ^V"ixQ)ix - xo)^ + ^V^^\x ~ xof + ... (2.99) 

where x is the particle coordinate. To zeroth order in the expansion, the ground state energy 
of the particle is Vq and its position is xq. The first correction to the energy is given by the 
inclusion of the quadratic term of Ea. H2.99|) and, when expressed in units of Vq, it scales as 

V"^^yVo - l/{Vo^^xo) where Xq is assumed to be also a typical length scale of V{x). The 
first correction to the mean position of the particle is given by the inclusion of the cubic term 
in Eg. H2. 99(1 as a perturbation to the harmonic potential and, when expressed in units of xq, 
scales also as 1/{Vq^'^xo)1 

Coming back to the quantum gas problem, one has to produce a third order expansion of 
the model Hamiltonian and treat the resulting cubic Hamiltonian as a perturbation to H2. 
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The derivation of is detailed in an appendix of [28], we give here only the result: 



90 
2 

2 



^^^^^ /^mP ^ ' 



' 1/2 
1/2 

. Po 



Sp 



1/2 

Po 

1/2 

Po,+j 



Sp+, {9+, - 9) 



h 

8to 



ESp f -1/2 . 1/2 1/2 . - 

— [Po ^Po -Po ^Po 
Po ^ 



1/2 



16r7i 



5P' . ^ 5p' 



5/2 

Po 



3/2 

Po 



A 



(2.100) 

6p 



There are then two main ways to calculate the correction to the mean density due to . The 
first one relies on finite temperature perturbation theory: to calculate expectation values in the 
density operator exp[— /9(i/2 + ^3)] to first order in H^^ one uses the imaginary time version 
of the time dependent first order perturbation theory: 



-P{H2 + H3) 



-PH2 



(2.101) 



One is then back to the calculation of expectation values of operators in the thermal state 
corresponding to H2, and Wick's theorem can be applied. 

The second way to include the effect of perturbatively is here simpler and was followed 
in [28] . It consists in writing the equation of motion of 6 in Heisenberg representation for the 
Hamiltonian H2 + -^3: with respect to the evolution governed by i?2, 'new' terms appear. One 
then takes the expectation value of this equation in the desired thermal state. The expectation 
value of the 'new' terms can be taken in the unperturbed thermal state, as they originate from 
the perturbation H^. Furthermore the expectation value of dtO vanishes in steady state, as 
proved in an appendix of [28], so that 







2m 



U + 350P0 - P 
gopl^^AB^ B + 



((5/5)3 - {B^B)2 



1/2 
Po 



B 



'%-2{P')2^dN^,^o 



(2.102) 



where the thermal average (. . .)2 is taken with the unperturbed Hamiltonian H2 and (. . .)3 is 
taken with the perturbed Hamiltonian H2 + to first order in H-^. The fact that po is a 
minimum of Hq imposes that the differential operator in the first line of Eq. (|2.1U2|I is strictly 
positive [1,24] so that this equation determines {5p)^ in a unique way. 



2.3 Applications of the formalism 

We present here without derivation simple applications of the expansion of the Hamiltonian 
performed in the previous subsection. For simplicity, we restrict to the spatially homogeneous 
case of a gas with periodic boundary conditions in a box of size L along each direction of space 
and we take the thermodynamic limit L +00 for a fixed value of the chemical potential. The 
derivation of the formulas given below and their extension to the case of a finite size, spatially 
inhomogeneous system can be found in [28]. 

2.3.1 Equation of state 

Thanks to Ea. (|2.1()2|l we can relate the mean density to the chemical potential to first order 
beyond the pure quasi-condensate approximation. As derived in [28] 

-^P+ I [iUk + Vufuk + Vk{Uk + 14)] (2.103) 
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where Uk = l/[exp(/3efe) — 1] is the mean occupation number of the Bogoliubov mode of wavevec- 
tor k and Uk,Vk are the amphtudes of the BogoUubov mode functions Wklr) = ?7feexp(ik • 
r)/L''/^ and Vk(r) = Vk exp(ik • r)/L''/^: 



Uk + Vk 



The corresponding BogoUubov eigenenergy is 



1/4 



2ra \ 2m 



2fi 



1/2 



(2.104) 



(2.105) 



In what follows we shall need the large k dependence of the zero temperature value of the 
integrand in Eq. H2.103|) : 

Vk{Uk + Vk)^-^. (2.106) 



Ea. H2.103|) is not totally satisfactory yet as it depends on the grid spacing I of our lattice 
model, both through the integration domain D defined after Ea. H2.42() and through the l- 
dependence of the coupling constant Qq. Let us check, for each value of the dimension d, that 
the I dependence disappears in the limit I ^, X. 

In ID this is clearly the case as both and the integral in Ea. (|2.1U3|) have a finite limit 
when I — > 0. One can then replace go by the ID coupling constant gi, and by ] — oo, +oo[. 

In 2D the integral in Eg. (|2. 11)31) diverges logarithmically in the limit I 0, because of 
the T = contribution to the integrand. We calculate the low-Z behavior of the integral by 
splitting the integration domain D in a disk of radius tt/1, over which the integral is performed 
in polar coordinates, and a complementary domain, over which the integrand is approximated 
by its asymptotic behavior Ea. (|2.106|l . We obtain 



In 



- 1 + 2 In 2 - 4G/7r + 0(ZV^^ 



(2.107) 



where G is Catalan's constant. Remarkably this compensates the logarithmic I dependence in 
the value of l/goj see Ea. (l2.61|l . We arrive at the equation of state of the 2D gas: 



P 



mil 



In 



dk 

(2^ 



{Uk + VkYuk 



(2.108) 



where 7 is Eulcr's constant and 02 is the 2D scattering length. This relation is identical to 
the result (20.45) obtained by the functional integral method in [18]. It allows to show that 
the necessary condition ^ 1 for our treatment to apply imposes hi{\/ pa\) ^ A-k at zero 
temperature. 

In 3D, the integral in Ea. H2.103|) also diverges. Using the identity Ea. H2.58() satisfied by the 
bare coupling constant g^ amounts to subtracting the high k behavior of the integrand, so that 
the I — > limit may be taken to give 



fJ- = P9 + g 



(2^ 



{Uk + Vkfnk + Vk{Uk + Vk) 



mfi 
Wk^ 



(2.109) 



In 3D our result coincides with the one of the usual Bogoliubov approach for a condensate. 
What might be surprising at first sight is that our results in 2D and in ID also coincide with 
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the ones predicted by a blind application of the usual Bogoliubov method, even if there is 
no true condensate! The same conclusion applies for the ground state energy of a gas of N 
particles, as shown in [28], and as could be expected from the equivalence of H2 with the usual 
Bogoliubov Hamiltonian. This fact that the usual Bogoliubov approach gives the correct result 
was commonly used in the literature in ID, see e.g. [5], and was justified by an asymptotic 
treatment of the solution based on the Bethe ansatz in [7]. With our extended Bogoliubov 
method, we reach this conclusion in a simpler, more transparent and more general way. 

2.3.2 Density fluctuations 

Our treatment relies on the assumption of weak density fluctuations. We then have to check 
that the relative variance of the density fluctuations is weak: 

4 = << 1. (2.110) 

Po 

Following the discussion of ^2.1.21 we expect this variance to be the sum of two contributions, 
one coming from the fact that Sp^ is not a normally ordered operator, and a second one 
involving the second order correlation function of the field, 

.g2(r) = (^t(i.)^t(o)^(o)v3(r)). (2.111) 

Explicit calculations are performed in [28]. In the thermodynamical limit one finds indeed 
that 

where : : is the standard notation to represent normal order (with all the on the left and all 
the "0 on the right) and 

(:(5p2(0) :)2 _ 2 f d'^k ..^ , ^^,2 



pI po Jd (27r)'^ 



[{Uk + VkYnk + VkiUk + Vk)]. (2.113) 



To ensure the condition ei <C 1, we require that pol"^ ^ 1 and that the normal-ordered 
contribution is small. The zero temperature contribution to Ea. (|2.113|l is found in [28] to be 
always smaller than 1/pol'^ as soon as Z < ^, see the table below in the column 'quantum term'. 
The same conclusion holds for the thermal contribution at temperatures ksT < p, cf. table. 
At temperatures ksT > p, the estimates produced in [28] are summarized in the table below. 
In 3D, the normal ordered contribution is automatically smaller than l/pol'^ since I < X. In 
2D and ID this is not necessarily the case; in practice it is convenient to adjust the value of I 
so that the normal ordered contribution is on the order of 1/pol'^. 

What is the link with the function 52 ? A careful calculation of 52 has to involve the 
correction to the mean density due to and leads to 

g^ir) ^p^ + 2p f 7^[(t/fc + Vkfuk + Vk{Uu + 14)] cos(k • r). (2.114) 

The difference between Ea. (|2.113|l and g2{0)/p^ — 1 therefore involves only the small difference 
between the pure quasi-condensate density po and the corrected mean density. As a conse- 
quence, the estimates at ksT > p also apply for the contribution of the thermal part to 52(0) 
and give conditions to have weak density fluctuations; in ID, we recover the condition obtained 
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in the classical field model, x ^ 1, and the asymptotic behaviour of (72(0), see Eg. (|1 .8411 . The 
quantum term leads to a divergence of 32(0) in 2D and 3D in the mathematical limit / — > 0, 
but this divergence is spurious: one should not forget that our treatment is applicable only if 
PqI"^ ^ 1. Finally, we note that the condition to have weak phase changes between two neigh- 
bouring cells of the grid is automatically satisfied when ei ^ 1: one has indeed the estimate 
ei ~ 1/pol'' [28]. 





quantum term 


thermal term, ksT < fi 


thermal term, ksT > /i 


d = 1 

d = 3 


1 

PoC 
ln(C/0 


(ksT/pr 


Po(, 

{kBT/p.)\n{kBT/p.) 




PoF , 


Poe 
1 




Poi^ 


7>o>^ 



Estimates for (: 5/5^(0) :)2/po- 



2.3.3 Coherence properties 

The last application of the formalism deals with the first order correlation function of the field: 

5i(r) = (Vi^(r)^(O)) ~ (pi/2(r)e^[^(°)-«»lpi/2(o)). (2.115) 

As detailed in [28] one first expands the two factors p^/^ up to second order in Sp. Then one 
calculates the thermal average with respect to the quadratic Hamiltonian H2, proving thanks 
to Wick's theorem identities like 

(e»Ae)2 = e-((A«?)^>2/2 (2.116) 

where = 6{0) — 9{r). Then one includes to first order corrections due to H3 using Ea. (|2.1(Jl|l . 
which has the only effect of replacing po by the mean density p. Finally one takes the continuous 
space limit I 0. This leads to the result: 

In Mr)/p] = ^J ^ [{Ui + V^W + (1 - cosk • r). (2.117) 

Remarkably this expression for the coherence function is easily related to the one predicted by 
the usual Bogoliubov theory for a Bose condensate: 



gi(r) = pexp 



5i no 



(2.118) 



In 3D, when a condensate is present, the argument of the exponential is small, so that the 
exponential may be expanded to first order: one recovers g\ ~ gf"^. In the infinite ID gas, 
or the infinite 2D gas at finite temperature, where no condensate is present, the Bogoliubov 
prediction (7f°^(r) diverges to —00 at large distances, whereas the prediction of the extended 
Bogoliubov theory tends to zero. 

The ID case is treated in some details in [28]. At zero temperature, one finds that g\ decays 
as a power law at large distances: 

5i(r)^p(-ij (2.119) 
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where the length scale is ri = e^~^^/4, 7 = 0.577 21 . . . being Euler's constant. This reproduces 
a result of [35] obtained by the path integral formalism. At finite temperature, the coherence 
function decays as a power law: 

ln[gi(r)/p] = ^ + K + o(l/r°°) (2.120) 

where o(l/r°°) is a function that tends to zero faster than any power law at infinity and the 
expression of the temperature dependent constant K is given in [28]. The coherence length 
coincides with the prediction of [18]; it also coincides with the classical field prediction H1.82|l . 
probably because the asymptotic behaviour of gi is controlled by modes of arbitrarily low 
wavevectors, and therefore of arbitrarily low Bogoliubov energy e^, for which the classicity 
condition ksT ^ Cfe is satisfied. 

An existing trap in the literature (see e.g. [36]) is to calculate separately the large r be- 
haviour of the quantum contribution to gi (that is at T = 0) and of the thermal contribution to 
gi (that is the only involving rifc): one find that the quantum bit behaves as a power law, and 
that to leading order the thermal bit behaves as an exponential, so that one may be tempted 
to conclude that the full gi behaves as an exponential times a power law, in contradiction to 
H2.120|) . What happens in reality is that the asymptotic expansion of the thermal bit involves 
a subleading power law contribution that exactly compensates the one of the quantum bit. 

One may wonder at which distances the asymptotic behaviour Eg. (|2. 120(1 is reached [37]. 
The result at a temperature k^T < /i is that, at this critical distance, the zero temperature 
asymptotics Ea. H2.119|) and the finite temperature ones Ea. H2.120|) should approximately give 
coinciding values, as justified in [38]. The critical length is therefore 

/crit - AV(27r20 (2.121) 

within logarithmic accuracy, and is actually the length over which a cross-over takes place, at 
a given temperature, from a T ~ power law behaviour of gi to a finite temperature fc^T < fi 
exponential behaviour of gi. 

As a final point, we emphasize that the present extended Bogoliubov approach gives access 
to the correlation functions of the field, like gi and 172, not only in the large distance regime, 
but also at a length scale on the order of the healing length ^. E.g. in ID it predicts that the 
third order derivative of gi in r = 0"*" is non zero, irrespective of the temperature: 

g^{^\0+) ^ fi^m^/{2h'^) (2.122) 

in agreement in the weakly interacting regime with an exact calculation at zero temperature 
based on the Bethe ansatz [39]. This is to be contrasted with other techniques like quantum 
hydrodynamics or Luttinger liquids [40,41], where only the large r part of the correlation func- 
tions is obtained, with the advantage however of not being restricted to the weakly interacting 
regime. 

3 INCURSION IN THE STRONGLY INTERACTING REGIME IN ID 

The first two sections of this lecture have presented simple methods to study weakly inter- 
acting degenerate Bose gases in arbitrary spatial dimension. In this last section, we present 
simple results on the opposite regime of extremely strong repulsive interactions, the regime of 
impenetrable bosons in ID, the so-called Tonks- Girardeau gas. The corresponding many-body 
problem can be studied exactly by a mapping on a gas of non-interacting fermions, a peculiarity 
of the ID case. 
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3.1 Mapping of impenetrable bosons onto an ideal Fermi gas 



This mapping is readily seen in the first quantization formaUsm for a gas of bosons with 
continuous spatial coordinates and interacting with a delta potential, gi6(xi — X2), in the 
limit gi — > +00 [42]. Such an infinitely repulsive interaction indeed imposes that the A'^-body 
wavefunction vanishes when two particles are at the same location. On the fundamental domain 
of coordinates of particles 1, . . . , TV sorted by ascending order, xi < . . . < xn, the wavefunction 
of a bosonic TV-body energy eigenstate is then an eigenstate of the kinetic energy plus trapping 
potential energy operator which vanishes at the border of the domain, and can therefore be 
shown to coincide with an eigenstate of A'^ polarized non-interacting fermions, which also 
vanishes when two particles are at the same location for the physically totally different reason 
that it is totally antisymmetric with respect to any permutation of particles [5, 7,43]. Out of 
the fundamental domain, the bosonic wavefunction and the fermionic wavefunction can differ 
by a global sign. As a consequence, the observables which involve the modulus squared of the 
A^-body wavefunction, like the pair correlation function (72, coincide for the bosonic and the 
fermionic wavefunctions; the observables that are sensitive to the phase of the wavefunction, 
like the coherence function gi or the momentum distribution, in general do not coincide for 
the bosonic and the fermionic wavefunctions. 

To perform explicit calculations, a second quantized version of the boson to fermion map- 
ping, the so-called Jordan- Wigner transformation, is actually convenient and easy to construct 
on a lattice model [40,44,45]. We therefore adapt the lattice model Ea. (|2.23|l to the case of 
impenetrable bosons in ID. Since go = +00, configurations where two bosons or more occupy 
the same lattice site are energetically suppressed, which amounts to introducing an orthogonal 
projector V on all the configurations with at most one boson per lattice site. At this stage, it 
is more convenient to manipulate the annihilation operator ax of one boson in the lattice site 
X than the field operator ip{x), the two operators being related by 



(3.1) 



where / is the grid spacing. The infinite go limit of the lattice Hamiltonian then reads 



E 



-—ValVMxV + {U{x) - ^i)ValVaxV 
2 m 



(3.2) 



In this form, the Hamiltonian is not easy to handle as Va^V and Val.V do not satisfy the 
usual bosonic commutation relations. In particular one has 







(3.3) 



{VaxV,ValV} = 1 



(3.4) 



which is reminiscent of a fermionic system. However it is not quite a fermionic system because 
the VaVs at different lattice sites commute rather than anticommute. One then introduces 



exp 



y<x 



Va,rV 



(3.5) 



or equivalently 



VcixV = exp 



iTT ^ clcy 



(3.6) 
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and one checks that the and cj, satisfy the usual fcrmionic anticommutation relations: 

{cx,Cy} = (3.7) 
{cx,cl} = 5x,y (3.8) 

After this transform, and if one uses the simplified representation Ea. H2.24|l of the discrete 
Laplacian, the model Hamiltonian Ea. (|3.2|) becomes 

= E -^^ct Ac, + {U{x) - ii)clc, (3.9) 

X 

which is a Hubbard model for non-interacting fermions. We used the fact that 

ValVax-iV = cle-'<-''-'cx-i = clcx-i- (3.10) 

One is left with a quadratic fermionic Hamiltonian which can be diagonalized [46]. In the 
subsections to come, we deduce some observables of the impenetrable Bose gas. 



3.2 Pair correlation function of impenetrable bosons 

The second order correlation function of the bosonic field, defined in Ea. (|2.111|l . is found to 
coincide with the one of the non-interacting fermions, since Val.axV — c\cx'. 

92{x) = ^{clclcocx). (3.11) 

The gas is at thermal equilibrium with temperature T. Wick's theorem can be used for the 
fermions to express 52 in terms of the first order coherence function f{x) of the fermions, 
and then in terms of the occupation numbers of the fermionic eigenmodes. For a spatially 
homogeneous system of density p and in the thermodynamic limit pL +00 and continuous 
limit pi ^ Q one finds 

i?2(x)=p2-|/(x)|2 (3.12) 



with 

27r exp[/3(n'^/cV(2m) - p)] + l 



+ 00 ^-^kx 



fi^)-l t;Z r...2,2... ^ (3-13) 



As expected for impenetrable bosons or for polarized fermions, .92(0) = 0. At zero temperature, 

/(.) = ^ (3.14) 

where kp = irp is the Fermi wavevector, so that the size of the hole of (72 is the mean interparticle 
separation \/p. In the non-degenerate regime T > Tp the hole size scales as the thermal de 
Broglie wavelength. 

An interesting application of 32 is to calculate the fluctuations of the number Nx of im- 
penetrable bosons in some sub-interval of length X of the bulk gas. The expectation value of 
Nx is pX and its variance is 

/• 

YaiNx = pX + dx dx' [g2{x - x') - p"^]. (3.15) 
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A more tractable formula involving only a single spatial integral is derived as follows: 

dx / dx'\f{x - x)\^ = / dxlx dx\f{x)f (3.16) 



and then one integrates by part in the integral over x, taking the derivative of /j^ ^ dx' and 
integrating the factor 1. This leads to 

l-x 

Ya.rNx^pX-2 dx {X - x)\f{x)f . (3.17) 
Jo 

At T = 0, an explicit expression is then obtained using the special functions Si(a;) and 
Ci(x): 

YaiNx = — \ [-TTkpX + cos{2kFX) - 1 + 2kFX Si(2fcF^) - 7 - \n{2kFX) + Ci{2kFX)] 

(3.18) 

from which the large kFX behaviour follows: 

where 7 = 0.57721 ... is Euler's constant. This results in extremely weak relative fluctuations 
of the number of particles, well below the Poisson limit for kFX ^ 1, since the dependence 
with X is logarithmic. 

At finite temperature, the Fermi distribution is a smooth, C°° function of the momentum 
and decreases faster than an exponential at infinity. One can then show, by repeated integration 
by parts in Ea. (|3.13|) . integrating exp(— ifca;) with respect to fc, that f{x) decreases faster than 
any power law at x = cxd. As a consequence, in the limit of a large interval, 



p- / dx\f{xf 



X+ / dx\x\\f{x)\^ +o{llX'^) (3.20) 



where the remainder o(l/X°°) decreases faster than any power law. From the Parseval- 
Plancherel identity one finds the more illuminating form 

Y&vNx - X / —n{k) [1 - n{k)] (3.21) 
J-oc 27r 

where n{k) is the Fermi distribution function of the ideal Fermi gas. This coincides with the 
usual grand canonical result for a system of total length X treated in the thermodynamic limit. 
In the non-degenerate regime, n{k) can be neglected as compared to 1 and the variance of Nx 
is the one of Poisson statistics, as expected for independent classical particles. 
What happens in the opposite regime T ^ Tp ? From the identity 

df,n{k) ^ Pn{k) [1 - n{k)] (3.22) 

and the low temperature expansion of the density of a ID ideal Fermi gas at fixed chemical 
potential [2]: 



1/2 



_ 1 /2m^\ 

^ TT 



(3.23) 
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one obtains the low temperature expansion 



YaiNx pX 



(3.24) 



where ksTp = h^kp/2m is the Fermi energy [48]. The fluctuations of the number of impenetra- 
ble bosons in a large interval is therefore subpoissonian at < T ^ Tp, but still proportional 
to X. We have plotted in figureElthe ratio of the variance of Nx and of its mean value, in the 
large X limit, as function oiT/Tp- 




T/T^ 



Fig. 5. In the limit of a large interval of length X, ratio of the variance and the mean value of the 
number of impenetrable bosons in this interval, or equivalently of the number of fermions for an ideal 
Fermi gas, as function of the temperature in units of the Fermi temperature. 

For a fixed temperature T ^ Tp, one can find the critical length Xc below which the 
variance of the number of particles has the zero temperature behaviour Ea. (|3.19|l and above 
which it acquires the finite temperature behaviour Ea. (|3.21|l . A first way of finding Xc is to 
use Eq.(20) of [40] giving an approximation to the function f{x) at T <C Tp: 



T2 



sin kpx 



4r/ sinh^(7rTfci.x/2TF) 



(3.25) 



One sees that f{x) differs weakly from its zero temperature value as long as TTTkp\x\ ^ 2Tp. 
One the contrary, if TrTkp\x\ ^ 2Tp, |/(a;)|^ is exponentially small so that one is left with 
Eq. H3.21|l . This gives a cross-over length 



kpXc 



2Tf 



(3.26) 



A second way of finding X^ is to equate Ea. (|3.19|) and Ea. (|3.21|) . Assuming ln2kpXc ^ 1, 
this gives the same order of magnitude for Xc as in the previous equation. A third way is to 
imagine that one can introduce a fictitious box of size X with periodic boundary conditions, 
containing a fictitious Fermi gas with the same chemical potential as the bulk system. This 
box gives wrong predictions at zero temperature, as the variance of Nx would then be exactly 
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zero, but gives indeed the correct result Ea. H3.21|) for finite T, when the thermodynamic Umit 
approximation is appUcable for the fictitious box, which occurs when 

fcsT > A£; (3.27) 

where AE is the energy separation between two consecutive energy levels in the box around 
the Fermi energy. One finds AE ^ 2TTTL^kp /mX which leads to the same Xc as before, within 
a numerical factor. 

3.3 First order coherence function of impenetrable bosons 

The first order coherence function in the general case of a spatially inhomogeneous system is 
defined as 

v) = {i^\y)i^{x)) = y(ata,). (3.28) 

As the expectation value is taken over a density operator tr such that VaV — cr, one can add 
factors equal to V in the above expression. This leads to an expression for gi in terms of the 
fcrmionic operators: 

g,{x,v) = -l-\clc,e'^^l-^''''^) (3.29) 

where we have assumed that x < y [49]. 

We now go through a sequence of transformations of Eq. I|3.29(l . First we rewrite it as 

, Tr [ctc^e^e^l 

91{^.V) = -1-' Vrn (^-^^^ 

where A and B are operators that are quadratic in the fermionic variables. B is equal to 
—jSHoo, where (3 — X/ksT and the Hamiltonian Hoa is given by Ea. H3.9|l . The operator A is 
ITT times the operator counting the number of fermions on the sites from x to y. We introduce 
the matrices in the lattice basis of the one-body operators corresponding to A and B: 



A = Y^Ac^ffcicp (3.31) 
B = ^S„^clc^ (3.32) 



Let us introduce a matrix C such that 



e^e^ = e^. (3.33) 
Then, as we show in the appendix ^ e^e^ = e'^ where the operator C is defined as 

C = Y,Ccf3cicp. (3.34) 

a, 13 

We also derive in this appendix the following expressions: 



Tr[ctc.e^] _ / 1 



(3.35) 



xy 



Tr [eC] V 1 + e"'^ 
Tr [e*^] = det (1 + e^). (3.36) 
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What remains to be done is to eliminate C in terms of A and B: this is easy since C appears in 
the final result only through its exponential, so that Eq.lj^23 can be used directly. A further 
simplification arises from the fact that A is proportional to a projector. The calculations are 
detailed in the appendix, we give the result: 

5i(a;,y)-4f-^) det(l-25) for x<y (3.37) 



Sap = {caCa) for X < a, /3 < y. (3.38) 



2^ VI -25,,^ 

where the matrix S is defined on the lattice sites in between x and y by the following thermal 
averages: 

Interestingly, for the spatially homogeneous case, in the continuous limit I 0, where matrices 
are replaced by operators acting on a functional space and the dispersion relation for fermions 
is quadratic, Ea. (|3.37|l becomes formally equivalent to the expression of gi given in [50] for 
impenetrable bosons in the Lieb-Liniger model. Finally, using the matrix identity M^^ — 
*comM/detM relating the inverse of a matrix M to the transpose of its comatrix, we obtain 
a writing equivalent to Ea. (|3.37|l : 

gi{x,y) = Yidet{2Sai3 - 5a,i3)\x<a<y,x<i3<y for x<y (3.39) 

that is in terms of the determinant of the matrix obtained by suppressing the first column and 
the last line of 25 — 1 . This coincides with [44] and, within a factor of two, with Eq. (8) of [40] . 
Extracting analytically physical information for gi is quite technical. The key result is the 
absence of long range order for the homogeneous system at zero temperature [51,52]: 

where A = 0.92418 ... is a numerical constant and kpX ^ 1. We refer to [53] for a recent 
overview of the various analytical results. An important consequence of Eq. (12301) is that the 
momentum distribution of impenetrable bosons at T = 0, which is the Fourier transform of gi, 
7r(fc) = cos{kx)gi{x), diverges as 2A[kp /{2Trk)]^^^ when fc — > 0. 

In figure El is plotted the result of a calculation of gi{X) based on a numerical evaluation 
of the determinant in Ea. H3.39|l at T = 0. We have actually plotted the product of gi{X) 
with (kpXY^^ to reveal that gi exhibits damped oscillations around the asymptotic formula 
Ea. H3.40|) : on a direct plot of the bare gi, these oscillations are not visible in practice. Note 
that these oscillations have a maximum roughly at the half integer values of pX and have a 
minimum roughly at the integer values of pX: this is in contradiction with the asymptotic 
series expansion of [52] but is reproduced by the analytical approach of [53]. 

Can we understand the emergence of such a l/X^^^ asymptotic behaviour of gi at zero 
temperature ? We present here an interpretation with no pretension of rigor. The first step is 
to realize that the power law of the long range behavior of gi is not changed if one eliminates 
the factor cj^Cx/l in the Ea. (|3.29|l . This is suggested by the calculations of [44], and we have 
checked this fact numerically. What remains to be understood is the long range behaviour of 
the function 

G{X) = (e"*^) = '^"S"" (3.41) 

n=0 

where Nx — X]^o ^z^^ operator giving the number of particles in the considered interval 

of length X, and 7r„ is the probability of having n particles in the interval. G{X) is then the 
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Fig. 6. Coherence function gi{x) of a zero temperature spatially homogeneous gas of ID impenetrable 
bosons in the thermodynamical limit, calculated numerically from Eg. (13. 391 . with pi = 0.005. g\ was 
multiplied by (kp^Y^'^ ■> where the Fermi wavevector kp is related to the density p by fc^ ~ vrp, to 
reveal that g\ oscillates around its asymptotic value. Dashed line: constant A of Eg. 13. 4011 . 



difference of the probability of having an even number of particles and the probability of having 
an odd number of particles in the interval. This can be estimated heuristically by assuming 
that the probability distribution of the number of fermions in the interval is roughly Gaussian 
for kpX ^ 1 [54]. One can then extend the summation for n to — oo since the standard 
deviation ax of the Gaussian, while being much larger than unity, is much smaller than the 
mean of the Gaussian. Using the Poisson formula 

-\-OQ +00 

E E (3.42) 



where /(/c) = / dxf{x)exp{—ikx) is the Fourier transform of an arbitrary function f{x), and 
restricting in the right hand side sum to the leading terms with n = and n = 1, one obtains 

G(X) ~ 2 cos(7rpX)e"''''"i/2. (3.43) 
The expression Eg. (|3. 19(1 of the standard deviation in the large X limit leads to 

GiX)^^^^cos{7rpX). (3.44) 

A numerical calculation confirms this result, with a coefficient 0.33 rather than 0.36. 

The conclusion is that the l/X^^"^ behaviour of the first order coherence function in the 
ground state of impenetrable bosons in ID reflects a property of the counting statistics of a 
ID zero temperature ideal Fermi gas in an interval of length X. 
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APPENDIX 



A WICK'S THEOREM 

Consider the following problem: for a density operator (here for bosons for simplicity, we shall 
come to the case of fermions later) 



P - 



Z 



exp 



where all z/q, are strictly positive, and the a^'s obey bosonic commutation relations: 

[da, d/s] = 0, 



calculate the expectation value 

/ = (Si . 



. 62,, ) = Tr 



hi. 



(A.l) 

(A.2) 
(A.3) 

(A.4) 



where the biS are arbitrary linear combinations of d^'s and dj^'s. We consider an even number 
of factors as the expectation value of the product of an odd number of 6's vanishes. 

A physical example is simply the ideal Bose gas at thermal equilibrium in the grand canon- 
ical ensemble, with — (3{ea — m)- Note that the density operator in the canonical ensemble 
does not reduce to (jA.ljl as it involves in addition the projector P/v on the subspace with a 
fixed total number of particles equal to N. 

The explicit calculation proceeds in four steps, in a derivation inspired by the lectures on 
statistical physical of Jacques des Cloizeaux at the Ecole norniale superieure of Paris, France, 
in 1988. 

First step: assume that &i is either d^Q or dj,^. Then there exists a number A such that 



pbi = Xbip. 



(A.5) 



One can check this identity using the Fock basis. E.g. one finds A = expz^ct^ for bi — aa„. By 
using the invariance of the trace by a cyclic permutation of the operators, one transfers hi to 
the extreme right and one then puts it through p thanks to (jA.5|l : 

(6i...62„) = A(62...S2„5i). (A.6) 

The following commutation relation allows to reintegrate the factor hi in the extreme left: 



2n 



[biM 



= ^ 61 . . .hj-i\bi,bj]bj+i . . .b2n- 
i=2 



(A.7) 



Since the commutator of two &'s is a pure number, we obtain 

2n 



(^1 



i=2 



A- 1 



'j+i ■ ■ 



(A., 
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Second step: Apply this general formula for two operators: 

(S,S,) = (A.9) 

This identity, which in particular implies the Bose formula, allows to rewrite the general formula 
in a simpler way, without any apparition of A: 

2n 

(61 . --hn) = Y,{hb,){b2 . ■■bj-ib, + i . ..b2n)- (A.IO) 

Third step: the linearity of (|A.10p with respect to bi implies that it holds even if bi is an 
arbitrary linear combination of da's and aj,'s! 

Fourth step: iterate the formula (|A.10p down to 2n = 2. So the most general expectation 
value / can be expressed in terms of expectation values of binary products, which are of course 
known from the Bose formula. 

The result is found to have a simple structure if one introduces the concept of a contraction. 
One collects the first factor bi with one of the other factors, that we call 6^: one says that one 
performs a contraction of bi with ba- One contracts the next factor available, bp, with one of 
the factors left, bj. Note that bp = b2 if a 2, otherwise j3 — 3. One repeats this process until 
all the factors are contracted. Then 

{bi...b2n)= {biK){hby)--.{hb^)- (A.ll) 

all contractions 

Let US count the number of possible ways of contracting the 2n factors. There are 2n — 1 
possibilities for the choice of the companion of bi , that is for the choice of b^ ■ There are 2n — 3 
possibilities for the choice of the companion of bp, that is for the choice of b^... Finally, there 
is a single possibility for the choice of the companion of 6^, that is for the choice of b^^. The 
total number of possible contractions is therefore 

(2n-l)x(2n-3)x...xl = ^. (A.12) 
Let us give an example for n — 2: 

{kkkk) = {bib2){kk) + {bik){b2k) + {kk){kk)- (A.13) 

What happens for fermions? The annihilation and creation operators now obey anticom- 
mutation relations 

{a„,4} = 5ap (A.14) 
{aa,dp} = 0. (A.15) 

As the number operator aj^ac is now bounded from above by unity, the coefficient Va is no 
longer restricted to positive values. The same technique as for bosons can be applied to derive 
the following recursive relation: 

(Si . . . kn) = Y,{'lY{kbj){k ■ ■ ■ ■ ■ ■ kn). (A.16) 

i=2 
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The factor (—1)-' appears because bj 'crossed' j — 2 factors. This gives rise to the Wick's rule: 

(61 . . . b2n) = ^ f^contractioniblba) {bpb^) . . . {b.^b^) . (A. 17) 

all contractions 

A sign ±, that we denote as Econtraction: is associated to each contraction. Each contraction 
actually defines a permutation of the 2n factors, performing the following mapping: 

(1,2,3,4,. .. ,2n-l,2n)^ 7,. (A.18) 

econtraction IS simply the signature of this permutation, that is (—1) to the power the number 
of transpositions of two elements required to realize the permutation. Let us give an example 
for n — 2: 

{bihkk) = {bib2)(hk) - (bik){b2k) + (bik)(b2h)- (A.19) 
B SOME USEFUL IDENTITIES FOR GAUSSIAISI OPERATORS 

Consider operators A, B and C that are quadratic in the fermionic field, in the sense defined 
by Eq. 3113. 3213. 34|l . The square matrices A, B and C of the associated quadratic forms are 
hnked by Ea. |X^ . 

B.l ABC identity 

We prove that the identity e^e^ — e*^ holds. The idea of the proof is to consider each 
exponential factor as an evolution operator and to show that the composition of two Gaussian 
evolution operators gives a Gaussian evolution operator. 
Let us introduce the formal time evolution: 

c„(t) = e^^c„e-^^. (B.l) 

The time derivative of Cq, (t) involves the commutator of B and Cq. that can be calculated from 
the fermionic anticommutation relations: 

■^Ca^ir) ^ -y Ba,[3CpiT). (B.2) 

dr ^ 

This differential system with the initial conditions Cct(O) — Ca is integrated formally in terms 
of the evolution matrix e"^"^. Specializing to r = 1 we obtain 



e^c„e-^ = 5](e-«)^^c^. (B.3) 



By using the same type of identity with B replaced with A or with C we obtain 

P.I 

7 7 

= e'^c^e-^ (B.6) 
where we used Ea. (|3.33|l . Similarly one can show that 



A S-t ~B -A C --i -C 

e e c' e e — e c' e 



(B.7) 
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Consider now the action of e on an Hartree-Fock state defined by an arbitrary list 
7i , . . . , 7„ of single particle states: 

|V')^ct^...ctjO). (B.8) 

We perform the following sequence of transformations: 

e^l^) = (e'^c;e-^)...(e^ct„e-^)|0) (B.9) 
= (e^e^c;e-^e-^)...(e^e^ct„e-^e-^)|0) (B.IO) 
= e^e^lV') (B.ll) 

where we used the fact that e-^\Q) = e--^|0) = e-'^lO) = |0) as can be checked by the power 
series expansion of the exponential. The operators e*^ and e^e^ have the same action on all 
the Hartree-Fock states, which form a basis, so they are equal. 



B.2 Trace and one-body expectation values of a Gaussian operator 

What is the trace of e"^ ? What is the one-body operator corresponding to e'-^ ? We give 
the answer assuming that the matrix C is diagonalizable, which is sufficient to conclude in the 
general case using a continuity argument. 

We call Ci the set of eigenvectors of C, with eigenvalue Ci, so that Cci = CiCt, and we call 
fj the adjoint basis normalized such that 

/; • e, = 5.,, (B.12) 

so that the closure relation is satisfied: 

i 

where {ei)a is the component of the vector in the original single particle basis. 
We introduce the corresponding creation and annihilation operators: 

Y.^o.{h)l (B.14) 
E^KeO- (B.15) 

Note that in general the creation operator 6^ is not the Hermitian conjugate of the annihilation 
operator hi since the matrix C is not necessarily Hermitian. This is why we introduced the 
symbol -I- rather than f in the exponent. However this has little effect on the algebra of 
second quantization! The hiS obviously anticommute among themselves, since they are linear 
combinations of the Cq's; also, their action on vacuum on the right gives zero. Similarly, the 
6^'s anticommute among themselves and their action on vacuum on the left gives zero. Finally, 
one derives from Ea. ljB.12|l the anticommutation relations 

{6„fe+} = <5,,,. (B.16) 

We can then construct a Fock basis as 

|M)^(6+)-...(6+)""|0) (B.17) 
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where M is the order of the matrix C, that is the dimension of the single particle state space, 
and the occupation numbers are equal to or 1. The adjoint Fock basis is 



(Mi^(oi(feMr^..(&ir- (B.18) 

The advantage of these algebraic constructions is that the operator C has now simple 
expressions and actions: 

C = ^C,6+&, (B.19) 

i 

and, from the usual calculus in second quantized form: 

C\{n})={^C.n}j\{n}). (B.20) 



As a consequence 



One has also 



Tr [e^] =^ne"'C. ^det + (B.21) 

{„} ^ 



Tr 



(B.22) 



Tr [e^^] 1 + e- 

and finally, using the closure relation Ea. (|B.13|l . we get Eg. 35(1 . 

B.3 Back to gi 

We assume now that the matrix A takes the specific form inQ where the matrix Q is a projector 
on the sites in between x and y: 

V 

Qa/3 = <5a,/3 ^ (B.23) 
z—x 

This leads to strong simplifications. First one has 

6-^ = 1- 2Q. (B.24) 

So, from Eq. 



1 + = 1 + (1 - 2Q)e^ = [1 - 2QS] (1 + e^) (B.25) 
where we have introduced 

5.^. (B.26) 

This last matrix has a very simple physical interpretation: using Ea. H3.35f) with C replaced by 
B, we find that Saf3 is the thermal average of cjjCa. We finally obtain 

T-^-e^T^ = (l-2Q)5 \ (B.27) 

The matrix to invert in the above expression can be put in block triangular form by collecting 
the lattice sites in between x and y in one block, and the other lattice sites in a second block. 
In block notations, one can check that 



blockii blocki2 \ " _ ( blockj^j^^ — block^^^ blocki2 block22 



block22 ) V block^2^ 



(B.28) 
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Here the right upper block of the inverse matrix is not useful as the lattice site y, on which the 
inverse matrix acts in the expression for gi{x,y), belongs to the first block. One has also that 
the determinant of a block triangular matrix is the product of the determinants of blockn and 
block22- This leads to 



When speciahzed to the case x ^ y, this gives Ea. H3.37|l . 
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